Recent change in spatial distribution of the European flat oyster (Ostrea edulis) inferred from field data and empirical models of living oysters and empty shells

Abstract Marine coastal areas are increasingly affected by human activities resulting in changes in species and habitat distributions. Understanding these patterns and its causes and consequences is important for conservation and restoration of such changing habitats. One habitat that has been heavily affected by human use are the North Sea oyster beds which once were abundant but have lost large parts of its coastal distribution due to overexploitation. Based on data of living and dead assemblages of Ostrea edulis collected using video transects, we used an ensemble modeling technique to model and predict current and recent distribution of O. edulis along the Swedish west coast where its distribution is, in relative terms, still rather unaffected. We could detect a recent change in the distribution of O. edulis along the coast which to a large extent could be attributed to a change in depth distribution, suggesting that the population of O. edulis have a slightly shallower distribution today than in the past. Although a potential mismatch between living and dead assemblages, caused by a complex combination of biological and environmental conditions, needs to be considered in the interpretations drawn, it may be a way around the lack of suitable background data in management decisions. This provides important information for management and conservation of the native oyster beds. Furthermore, this study illustrates a method for identifying recent changes in species distribution using dead assemblages of bivalves.


| INTRODUC TI ON
As human population have grown, so have the pressure on coastal environments. Humans have long been dependent on the coastal zones and the services its ecosystems provide (Mehvar et al., 2018;Neumann et al., 2015). Thus, human activity has highly modified marine ecosystems in temperate regions and, in combination with effects of global climate changes, many of these ecosystems are far from their natural and original status (Halpern et al., 2008(Halpern et al., , 2015, and the speed at which the pressures on the ocean increases are growing (Halpern et al., 2019). One of these intensively used areas is the North Sea which, in a relatively short time, has lost almost all its offshore and large part of its coastal oyster grounds due to overexploitation (Beck et al., 2011). Once abundant and ecologically important, reef habitats have vanished from the ecosystem and today, beds of the European flat oyster (Ostrea edulis) are rare or absent from most of their natural range (Airoldi & Beck, 2007;Haelters & Kerckhof, 2009;Laing et al., 2005). These reef structures are now protected and should, according to the Habitats Directive, be preserved and/or restored to favorable conservation status.
Oyster reefs produce a variety of important and valuable ecosystem services: improved water quality, increased nutrient uptake, complex three-dimensional structure which provide habitats, food, and protection for a large number of species, as well as providing a valuable food and economic resource for humans (Gerritsen et al., 1994;Grabowski et al., 2005;Nelson et al., 2004).
Restoration efforts of Ostrea edulis are currently restricted to the few remnant populations which is already heavily impacted by historical exploitation. Since the "natural" state of oyster beds have often been defined at the time when the protective area was designated, the baseline upon which changes are measured and on which management is determined is often already an affected state which might not reflect the past conditions. Assessing fully natural conditions is often difficult as impacts most often pre-date detailed documentation of the areas. Thus, conservation, restoration, and management are often based on incomplete knowledge and reference conditions. However, due to the properties of oyster beds inherited by the fact that the oysters produce relatively thick and stable shells which remain long after the animal itself dies, dead shell assemblages is considered one of the most crucial components of a healthy oyster reef and is thus considered a main component in oyster reef restoration efforts. Comparison of living and dead assemblages offers possibilities to study how processes of settlement and decomposition influence the preservation of species in fossil record (Kidwell, 2013), although caution is needed when interpreting results as anthropogenic impacts are known to cause mismatch between living and dead assemblages (Kidwell, 2013;Kidwell & Tomasovych, 2013). However, the dead shells also offer another less studied opportunity, the potential to back track recent changes in the spatial distribution of oyster beds. Applications of dead assemblages of organisms in investigations of ecological changes have been scarce, but these studies have demonstrated that the distribution of dead organisms can function as a long-term record of changes in community structure and function (Liversage et al., 2020). For example, Kidwell (2007Kidwell ( , 2008 investigated community changes in response to anthropogenic activities (e.g., fishing) using dead assemblages.
Using dead assemblages for tracking short-and long-term changes in spatial distribution of the living community over time can potentially offer an important opportunity to follow changes in ranges due to, among others, human pressure when time series of a species large-scale spatial distribution is lacking. How far back you can go depends on species and local biological and environmental processes. For oysters, this is limited by the destruction and dissolution of its shell which in turn depends on a combination of biological, geochemical, and sedimentary factors; the same processes that ultimately control the rate of loss from natural reefs. Estimates of halflife times for oyster shells vary from a few years to a few decades. Powell et al. (2006) calculated the half-life time for Eastern oysters (Crassostrea virginica) to between 1 and 20 years while Waldbusser, Steenson, et al. (2011) estimated it to be up to 40 years. The shells might also be transported away from its origin contributing to the uncertainty in the time span available. However, this transport is not well known and thus difficult to account for, and with rather heavy shells, of oysters, is restricted to more exposed areas. With shellboring organisms generally considered to be the primary reason for oyster shell degradation (Carver et al., 2010), longer half-life times are expected in areas scares of such organisms.
With humans increasingly altering the distribution of species by modifying habitat, changing global climate, and introducing new species (Mack & Lonsdale, 2001;Parmesan, 2006;Tilman & Lehman, 2001), successful management of biological resources, such as oyster or mussel beds, depends on our ability to predict the potential spatial distribution of range-changing species and to understand the forces that limit their distribution. Combining knowledge of distribution of living and dead oysters (or other bivalves) with species distribution modeling (SDM) and geographic information system (GIS) offers the possibility to predict and map current spatial distribution (based on living oysters) and former recent areas (based on dead shells), which can be used in conservation to identify areas of special interest for restoration and conservation project.
The aim of the study was to investigate similarities and differences in predictive models and predicted distribution of adult living Ostrea edulis densities >1 ind.m −2 and sites with high abundances (>1 ind.m −2 ) in dead shells to identify and evaluate recent changes in the spatial distribution of the population along the Swedish west coast. We also demonstrate a potential new method for predicting short-term changes in species distribution based on species distribution modeling technique and field data on living and dead shells.
For these purposes, we generated two models, using a similar approach as Bergström et al. (2021), to (a) predicted distribution of high-density areas of living oysters (LO), and (b) predicted distribution of areas with high densities of empty oyster shells (EOS) in the absence of high densities of living oysters. These were used to identify differences between areas with living high densities of O. edulis and areas in which the abundances of shells indicate a recent, but no longer, high abundance of living oysters, and finding places with either living or dead high densities. We also used the results to discuss potential explanations to changes in distribution and potential distinguishing characteristics of sites which previously have had high densities (i.e., those with a lot of dead oysters) but which now have none.

| Study area and species distribution data
This study was carried out in the northern part of the Swedish west coast bordering Skagerrak ( Figure 1). This area is dominated by a mixture of rocky and sandy shores spread across an extensive archipelago with large number of islands, of which the inner parts may be ice covered during the winter, and it host the vast majority of the Swedish native flat oyster population. A dataset of abundance of living Ostrea edulis and dead shells was extracted from the studies performed by Thorngren et al. (2017Thorngren et al. ( , 2019. Thorngren et al. (2019) visited a total of 452 randomly selected locations in the study area during a 2-year period (2013)(2014). Sampling was restricted to a maximum depth of 10 m and stratified into three depth strata (0-3, 3-6, and 6-10 m) and to areas classified as moderately exposed or less according to the classification of wave and wind exposure (Naturvårdsverket, 2006). In each of these locations, two 20-m-long, 0.8-m-wide video transects were filmed, at a speed of approximately 0.4 knots, using a downward-facing high-definition camera mounted 50 cm from the sea floor to a towed sledge. The video films were then analyzed according to Thorngren et al. (2017) with number of adult sized "living," "possibly living," and "dead" oyster recorded.
From the results presented in Thorngren et al. (2017), which investigated small-scale variability, we know that the distribution of living Ostrea edulis is patchy at the level of individual transects. But, having long enough transects, the effect of this patchiness is negligible in this study. This sampling method was validated in Thorngren et al. (2017) by comparing video data, from filmed transects, with field measurement made by snorkeling in the same transects. The overall correct classification rate using this method was found to be 0.81, while the significant variability between observers were less than 1%. Detailed results are presented in Table 2 and the associated text in Thorngren et al. (2017).
Although often a density of 5 ind. m −2 is conventionally used in conservation purposes, e.g., OSPAR, to define oyster beds (Haelters & Kerckhof, 2009), we used a more inclusive level of 1 ind. m −2 as definition of high-density areas in this study. The rationale behind this is the observation made by Thorngren et al. (2019) and Bergström et al. (2021), which showed that high-density areas using the OSPAR definition would correspond to roughly 1% of the visited sites while the one selected here would include roughly 5% of the sites and 85% of the oyster population. Thus, the selected level of 1 ind. m −2 provides a better separation between those sites that contribute most to the total population from the rest while at the same time provide a solid base for modeling. For further details on experimental design and sampling, see Thorngren et al. (2017Thorngren et al. ( , 2019. F I G U R E 1 Overview of study area (a) and sample intensity (b) with close-up on the most intensive sample areas (c) within the Kosterhavet National Park

| Environmental data and pre-modeling analysis
Three environmental predictor variables (depth, salinity, and exposure) were selected for the modeling. The candidate predictors were based on the general relevance for species distribution in marine environment, their availability as full covering raster layers for the investigated area, and the results obtained by Bergström et al. (2021). Wave exposure data from the Naturvårdsverket (2006) were corrected for depth according to Bekkby et al. (2008)  field. Pre-modeling, the selected explanatory variables were tested from correlation using variance inflation factor analysis. This is a simple approach to identify collinearity among predictor variables and those who obtain a high value in the calculations are removed.
Although arbitrary, a value of 5-10 is normally considered high collinearity (Akinwande et al., 2015;Hair et al., 1995;Kline, 1998), here we used a more conservative level of 3. These initial analyses showed that the tree environmental variables were largely uncorrelated (VIF value below 3) within the sampled data and were thus considered sufficiently uncorrelated for use in modeling.

| Modeling algorithms
We fitted an ensemble species distribution model based on nine commonly used algorithms: four machine learning methods, Artificial Neural Networks (ANN, Ripley, 1996), Classification Tree Analysis (CTA, Breiman et al., 1984), Generalized Boosted Models (GBM, Ridgeway, 1999), and Random Forest (RF, Breiman, 2001); four regression-based methods, Generalize Additive Models (GAM, Hastie & Tibshirani, 1990), Generalized Linear Models (GLM, McCullagh & Nelder, 1989), Multivariate Regression Splines (MARS, Friedman, 1991), and Flexible Discriminant Analysis (FDA, Hastie et al., 1994); and one envelope-style method, Surface Range Envelope (SRE, Busby, 1991). In this study, the BIOMOD2 package (Thuiller et al., 2016) for R software (R Core Team, 2019) and its related packages were used for all the selected algorithms. The default settings for modeling options in BIOMOD2 were used which do not include interactions for the regression-based methods. The reason for this was that incorporating interaction terms would quickly increase the number of effective variables. A consensus ensemble approach was applied using the BIOMOD2 platform models generated by the selected individual models. This approach should in theory offer a more robust, less noisy, prediction for the potential and realized distribution of Ostrea edulis than single algorithm models.
Models were built using a 100-fold cross-validation and randomly splitting the data into training (70%) and test data (30%) for respective model calibration and testing. This splitting procedure permit evaluation of model accuracy and predictive performance when data are non-independent.

| Modeling evaluation
The performance of the models (i.e., strength of agreement among distribution data and each model) was assessed using the area under receiver curve (AUC, Hanley & McNeil, 1982) and true skills statistics (TSS, Allouche et al., 2006;Liu et al., 2005) together with sensitivity (percentage of good presence predictions) and specificity (percentage of good absence predictions). AUC is threshold independent and evaluates both false-positive error rate and the true-positive rate in order to obtain a measurement of the model accuracy while TSS takes both omission and commission errors into account (Allouche et al., 2006) maximizing the sum of sensitivity and specificity. AUC range from 0 to 1 with values below 0.5 representing models that is not better than random and a value of 1, a highly accurate model (Scarnati et al., 2009). Although AUC has been highly criticized in some studies (Austin, 2007;Jimenez-Valverde, 2012;Raes et al., 2009), it is still the most commonly use measure to assess model accuracy and therefore considered a useful measure for this study. TSS range from −1 to 1 with values above 0.6 considered useful (Coetzee et al., 2009). To ensure an accurate model prediction, only individual models above a critical TSS (0.6) value were implemented in the final ensemble model, further the included models were weighted based on their TSS values. This was done because weighted means have been suggested to have the best performance of the ensemble methods available (Marmion et al., 2009).
Functional relationship between the explanatory variables and the response, either occurrence of high-density (>1 ind. m −2 ) living oysters (LO) or high-density (>1 ind. m −2 ) dead shells (EOS), was further explored using partial dependence plots. The predictors that most strongly influence the model were determined using variable importance assuming that the most important variables are the ones with a relative importance above the mean of the predictor variables within each subset. Variable importance was estimated using the built-in function in the BIOMOD2 package (Thuiller et al., 2016), which is based on shuffle a single variable of the given data, making model predictions using this "shuffled" data and then computing the Pearson correlation between reference and "shuffled" predictions returning a value between 0 and 1 were the higher value the higher importance of that variable. One limitation of this method is that it does not account for interactions between different variables.

| Prediction and post-modeling analysis
To evaluate the potential recent loss of high-density oyster habitats, after the final ensemble models were created and evaluated for variable importance and predictive power of the environmental variables, they were applied to high-resolution (15 × 15 m) raster layers of the environmental variables to generate probability maps of occurrences (living oysters and dead shells, respectively). These probability maps were translated into presence-absence distribution using the threshold (cut-off, optimized based on a data-driven approach using Youdens index) calculated by BIOMOD2 where all areas (15 × 15 m raster cells) with a predicted probability above the threshold (>cut-off) grouped into "present," whereas lower suitability values were grouped into "absent." The observed results were then analyzed visually to identify areas of high interest for further studies on changes in the distribution of living oysters along the Swedish coast. Since areas based on cut-offs optimized using Youdens index (minimizing the total error, that is, false negatives plus false positives) for each model are not easily comparable between models, we also used a more "decision-analytic approach" where the same cut-off values were used for both models. This allowed us to estimate areal extent of LO and EOS plus overlapping areas, under conditions where the criteria for classification are equal for LO and EOS models.

| Pre-modeling results
Of the total 452 sites investigated, 59 sites had high densities

| Model performance and variable importance
The  in the reverse order of importance (Table 2). However, in the model

F I G U R E 2 Correlation and variation inflation factor among the three environmental predictors in the model
including empty shells (EOS), the difference between the importance of depth and the other predictors is smaller and depth and salinity had almost the same importance. Although more important for EOS model, the partial plots (Figure 4), showing the functional relationship between environmental predictors and the probability of occurrence showed the same general pattern for both EOS and LO models with regard to salinity and exposure (Figure 4b,c,e,f).
However, with a slightly more distinct peak in partial dependence of exposure for EOS models than LO at the middle of the investigated exposure range.

| Prediction and post-modeling results
As illustrated by the partial plots in Figure 4 and the examples in Furthermore, although sites of only dead assemblages were found throughout the entire northern part of the study area, the highest frequency of sites with only dead assemblages was observed in two areas close to the northern and southern limits of the "high-density area," respectively. The placement of these areas with dead assemblages indicates that the total spatial distribution of high-density areas has shrunk slightly ( Figure 6). Figure 6 also illustrates that there seem to be slightly less occurrences of sites with only dead assemblages in more sheltered areas compared to more exposed areas.
In order to estimate the overall areal extent and map the locations of living oysters (LO) and empty shells (EOS), we utilized predictions of high densities of living oyster and empty shells. The areal extents were predicted using both model-specific optimized cut-offs (Youdens index) and the decision analytic approach of equal cut-offs. These analyses showed that the predicted areal extent of living oysters and empty shells was largely sensitive to cut-off values (Table 3, Figure 5). Using model-specific optimization to minimize total error, we obtained an overlap of only 3% between living oysters and empty shells of the total predicted area while living oysters covered 75% of the area. Instead, if using equal sized cut-offs, the overlap was 10% and the dead assemblages made up 45% of the total predicted area ( Table 3). The corresponding observed values were 34% overlap, 12% areas with only living assemblage, and 54% of the total observed sites with high densities having only high densities of empty shells.

| DISCUSS ION
For meaningful targets for conservation and restoration, it is vital with baseline information on the status of the community/species before a recent and/or rapid change (NRC, 2005  For environmental assessment, the mismatch between living and dead assemblages is not perfect, as there are some overlap between dead assemblages and pristine areas and the assemblages of dead shells may represent a range of time spans depending on local biological and environmental conditions. Often these assemblages of dead shells include a majority of relatively recent deaths and a F I G U R E 5 Examples of predicted distribution of high densities of living oyster and empty shells showing a general pattern of changed distribution of living oysters. The pairs of maps (a-d, b-e, and c-f) represent the predictions using optimal cut-offs (a, b, and c using Youdens index) and a more "decision-analytic approach" using equal cut-offs (d, e, and f). Land is presented in green, Predicted areas of high densities of living oysters in blue and high densities of dead assemblages in striped long tail of rare shells from individuals that died at increasingly earlier time. Given these time spans, the distribution of dead shells may not be a perfect indicator of present or past habitats, restricting the interpretations and conclusions that can be drawn from the mismatch in distribution between living and dead assemblages (Powell et al., 2017). Additionally, potential accumulation of transported shells further complicates the interpretation (Callender et al., 1992;Miller, 1988;Zenetos, 1990) as will the tendency of species not to be distributed in all suitable habitats during all time (Levinton, 1970;Powell et al., 1986). Thus, the spatial distribution of dead shells and living animals may not always be a suitable indicator of present and/ or past habitat. However, given the lack of background data available in most areas and the urgent need for baseline data to support management and conservation decisions, the use of living versus dead assemblage's might potentially be the only way to provide evidence of recent changes in distributional patterns in the absence of widespread baseline information. Tomašových and Kidwell (2009) showed that death assemblages have the same ability to capture environmental gradients as living assemblages, further supporting the idea of using dead assemblages to track recent changes in species distributions. This potential is to a large degree unaddressed in marine environments. Furthermore, in cases when historical and contemporary knowledge on the vertical distribution of regional species populations is inadequate, as for the Swedish Ostrea edulis population, the use of methods utilizing dead assemblages might be the only way to obtain some scientifically based baseline information for future management of the population. As in the case of Ostrea edulis, heavier shells require a much higher water movement energy for transportation, making them more suitable for utilization of recent changes in spatial distribution than lighter shells such as blue mussel (Mytilus edulis) shells which are more easily transported away from its origin. The observed increased importance of salinity and exposure for models of the distribution of dead shells is not fully matched by an alteration in partial dependence suggesting that increased impor-

TA B L E 3
Comparison of occurrence of high-density areas with living oysters (LO), empty oyster shells (EOS), and both LO and EOS for observed and modeled data of medium exposed sites for the living assemblages. One potential explanation to this increased importance could be that the dissolution rates of Ostrea edulis shells are altered by salinity and exposure. Another explanation would be that shells of oysters growing in lower salinities are thinner due to higher cost of calcification (Malone & Dodd, 1967;Waldbusser, Voigt, et al., 2011) and thus dissolved faster upon death of the oyster. This increased cost is partly a result of lower salinity causing reduced Ca 2+ concentration and lower levels of total inorganic carbon in the water (Hofmann et al., 2009;Mook & Koene, 1975). However, this is a less likely explanation, as this would have resulted in changed patterns in the partial dependence plots between the EOS and LO models, which is not observed (Figure 4). As the oyster shell is a dynamic resource subject to a number of processes causing degradation, it is most likely a combination of several processes explaining this observed difference and untangling these are beyond the scope of this study.
One candidate process is increased water turbidity and local resuspension which are known to affect the growth of O. edulis with decreased growth, as the beneficial effects of resuspension (e.g., increased food supply) are inhibited by increased resuspension due to decreased ingestion and dilution of food with inorganics (Grant et al., 1990) and inhibiting recruitment (Moore, 1977).
The reason behind the observed pattern ( Figure 6) of high frequency of dead assemblages in mainly two areas is currently unknown and needs further studies. Hypothesis about the causes includes loss of genetic variability inflicting weaker resilience toward infrequently occurring environmental events and diseases.
Another would be local outbreaks of diseases or parasitic events affecting local populations. One interesting note is that the two areas Species distribution models are an increasingly important tool in conservation and management (Hao et al., 2019). Using several models combined in an ensemble modeling approach may be a safer and more reliable method that can overcome uncertainty in model selection, thus further improving the reliability of the predicted spatial distribution.
Using species distribution models to further evaluate the spatial distribution of living and dead assemblages might be a possible way to estimate recent change in potential spatial distribution for larger areas. However, its interpretations require a clear distinction between potential and realized distributions (see Soberón, 2007).
Using, as in this study, field data on living and dead shells will give models on the potential distribution (i.e., places where the spe-  Figure 5). Because the optimal cut-off for living and dead oysters (based on a data-driven approach using Youdens index) differs strongly (LO ≈ 0.1 and EOS ≈ 0.3), areas were potential not easily comparable between the two models. While it might be argued that the main concern in this context of this study is error due to false negative, we also used a more "decision-analytic approach" (Steyerberg et al., 2011) where cut-offs values were set equal for both models. This allowed estimation of areal extent of living and dead assemblages plus overlapping areas under conditions where the criteria for classification were equal. Therefore, due to the sensitivity of predicted areal extent to cut-off values and the approaches used, caution is needed when evaluating and drawing conclusions from comparisons of models, and the conclusions might change depending on either data-driven approach or more decision-analytic approach is used.
Prediction based on models relies on full-covering spatial data of all the predictors used, which are not always available. Naturally, several other environmental parameters such as pH, temperature, food availability, and gravel content in sediment (Bayne, 2017;Bergström et al., 2021;Cano et al., 1997;Davis & Calabrese, 1969;Laing & Spencer, 2006;Laing et al., 2005;Pardo et al., 2016) have

CO N FLI C T O F I NTE R E S T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data from Thorngren et al. (2019)